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ABSTRACT 

We demonstrate that the observed increase of some nebular line ratios with height 
above the midplane in the diffuse ionized gas (DIG) in the Milky Way and other 
galaxies is a natural consequence of the progressive hardening of the radiation field 
far from the midplane ionizing sources. To obtain increasing temperatures and line 
ratios away from the midplane, our photoionization simulations of a multi-component 
interstellar medium do not require as much additional heating (over and above that 
from photoionization) as previous studies that employed one dimensional, spherically 
averaged models. 

Radiation leaking into the DIG from density bounded H ii regions is generally 
harder in the H-ionizing continuum and has its He-ionizing photons suppressed com- 
pared to the ionizing source of the H ii region. In line with other recent investigations, 
we find that such leaky H ii region models can provide elevated temperatures and line 
ratios, and a lower He"*" fraction in the DIG. For a composite model representing the 
relative spectral types of O stars in the solar neighbourhood, we find that additional 
heating less than 10~^^no ergs/s/cm'^ can reproduce the observed elevated line ratios 
in the DIG. This additional heating is considerably less than previous estimates due to 
the natural hardening of the radiation field reaching large heights in our simulations. 
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1 INTRODUCTION 

The presence of extended layers of diffuse ionized gas (DIG) 
in the Milky Way and other galaxies is inferred from faint 
emission in Ha and other nebular lines (e.g., Reynolds 1995; 
Rand 1998; Hoopes, Walterbos, & Greenawalt 1996; Dom- 
goergen & Dettmar 1997; Otte & Dettmar 1999; Hoopes & 
Walterbos 2003; Wang, Heckman, & Lehnert 1998; Otte et 
al. 2001, 2002). The ionization source for the DIG is be- 
lieved to be O and B stars close to the midplane of galaxies. 
Photoionization models using OB stars can reproduce the 
gross features of the DIG ionization structure (e.g.. Miller 
& Cox 1993; Dove & ShuU 1994). If extra heating is in- 
cluded, in addition to that from photoionization by O stars, 
ID models can reproduce the DIG spectrum (Mathis 1986; 
Domgoergen & Mathis 1994; Mathis 2000; Sembach et al. 
2000). The main problem with O stars as the ionization 
source of the DIG is that for a smooth density distribution 
the mean free path of Lyman continuum photons is very 
small, ~ 0.1 pc for n(H°) = Icm"^. Thus, it is difficult 



for photons to traverse the large distances to ionize the ex- 
traplanar gas which has a scaleheight of around 1 kpc in 
the Milky Way (Haffner, Reynolds, & Tufte 1999). How- 
ever, three dimensional photoionization models show that a 
two component (Wood & Loeb 2000) or fractal ISM (Ciardi, 
Bianchi, & Ferrara 2002) can provide low density paths al- 
lowing Lyman continuum photons to reach large \z\ heights 
above their midplane sources. 

Measurements of emission lines in addition to Ha pro- 
vide probes of the abundances, temperatures, and densities 
in the DIG. Some of the most studied lines in the Milky 
Way and other galaxies are [N II] A6584, [S II] A6717, [O HI] 
A5007, [O II] A3727 doublet, and [O I] A6300. Observations of 
the He I A5876 line provide information on the spectrum of 
the ionizing sources (Reynolds & Tufte 1995). Hereafter the 
above lines will simply be referred to as [N II], [S II], [O HI], 
[O II], [O I], and He I. Some differences in the observed line 
ratios in the DIG compared to traditional H ll regions are: 
[S II] /Ha and [N II] /Ha ratios increase with height above 
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the plane; [S II] /[N II] is quite uniform with height (Haffner 
et al. 1999) and with latitude (Rand 1997); He is observed to 
be underionized with respect to H (Reynolds & Tufte 1995; 
Heiles et al. 1996). 

Most models for the DIG employ smooth, one dimen- 
sional density distributions, and predict volume averages of 
the line strengths and ratios (e.g., Mathis 1986; Domgoergen 
& Mathis 1994; Collins & Rand 2001; Sembach et al. 2000; 
Mathis 2000). These models generally fail to reproduce the 
observed line ratios that increase with \z\ above the mid- 
plane. Haffner et al. (1999) and Reynolds, Haffner, & Tufte 
(1999) showed that the observed [S 11] /Ha and [N II] /Ha 
line ratios may be explained if the gas temperature increases 
with height above the midplanc. Including additional heat- 
ing over and above that of pure photoionization can reconcile 
the one dimensional models with observations (Reynolds et 
al. 1999; Mathis 2000). Additional heating may plausibly 
arise from photoelectric heating from dust (Reynolds & Cox 
1992), dissipation of turbulence (Slavin et al. 1993; Minter 
& Spangler 1997), and shocks (Raymond 1992). However, it 
is well known that in photoionized H II regions the radiation 
field hardens towards the edge of the Stromgrcn sphere re- 
sulting in the highest temperatures occurring at the largest 
distances from the ionizing source. This arises because low 
energy photons have relatively short mean free paths and axe 
absorbed close to the source. Higher energy photons have 
longer mean free paths, travel farther, and deposit more en- 
ergy per photon at large distances from the source, giving 
rise to the increasing temperature away from the source (Os- 
terbrock 1989). 

Such a scenario as described above is almost certainly 
occurring in the DIG, with the ionizing spectra penetrat- 
ing to large \z\ above the plane being significantly harder 
than the source spectra, naturally producing a temperature 
profile that increases with \z\. Models presented by Bland- 
Hawthorn, Freeman, & Quinn (1997), Wang et al. (1998), 
and Rand (1998) take this radiation transfer effect into ac- 
count by introducing a hardening of the radiation field with 
\z\. In these plane parallel models, the hardening of the ra- 
diation field leads to increasing temperatures at large dis- 
tances from the illuminated face and corresponding increases 
in [S II]/Ha and [N II]/Ha. Mathis (2000) also mentions this 
effect, stating that line ratios for lines of sight that pierce 
the outer edges of spherical models may reproduce the ob- 
servations. 

Radiation leaking into the ISM frtjni traditional H ll 
regions may also be hardened compared to the source spec- 
trum. Hoppes & Walterbos (2003) investigated photoioniza- 
tion by photons from leaky H ll regions finding that the 
hardened spectrum leads to elevated temperatures and in- 
creased line ratios when compared to models that do not 
envoke hardening of the source spectrum. In addition, leaky 
H II region models lead to a suppression of He-ionizing pho- 
tons {hf > 24.6 eV) and a corresponding decrease in ionizar 
tion stages such as He^ and N^"*". 

In this paper we present Monte Carlo photoionization 
models for a multi-component ISM. In what follows we de- 
scribe the photoionization code, adopted ISM density struc- 
ture, and spectra for the ionizing sources. Due to the very 
large parameter space, we restrict this paper to two dimen- 
sional models of a single source ionizing a multi component, 
stratified ISM. Note that although our models are for 2D 



systems, our simulations are run on 3D density grid. A fu- 
ture paper will present models for 3D geometries and illumi- 
nations, investigating the role of diffuse ionizing radiation, 
3D ionization and temperature structures, and the result- 
ing intensity maps. Some specific issues we address in this 
paper are: increased temperatures with increasing \z\ above 
the midplane; problems arising with fitting [S II] emission 
due to undetermined dielectronic recombination rates; pre- 
dictions of excess [O I] emission compared with observations; 
the ionization of He within the DIG. 



2 MODELS 

Our simulations use a 3D Monte Carlo photoionization code, 
with inputs being the ISM density structure, and the loca- 
tions and ionizing spectra of point sources within the simu- 
lation grid. These ingredients are briefly discussed below. 

2.1 Photoionization Models 

The photoionization simulations are performed using the 3D 
Monte Carlo code of Wood, Mathis, & Ercolano (2004). The 
code simulates photoionization due to multiple point or ex- 
tended sources and discretizes the density onto a 3D linear 
Cartesian grid. For all simulations in this paper the grid 
comprises 65"^ cells. The code performs well compared to tra- 
ditional ID codes and other independently developed Monte 
Caxlo codes (Ercolano et al. 2003; Och, Lucy, & Rosa 1998). 
We do not consider the effects of shocks or ionization fronts. 
In our simulations there is no C° or S" because we assume 
C and S are fully ionized by the ambient interstellar radia- 
tion field. There are no ionizing photons with energies above 
54.4 eV, so there is no He^"*", consistent with observations of 
almost all H ii regions. For further details of the code and 
comparisons with other codes see Wood et al. (2004). 

For most simulations, we adopt the following abun- 
dances by number relative to H: He/H = 0.1, C/H = 
140 ppm, N/H = 75 ppm, 0/H = 319 ppm, Ne/H = 
117 ppm, and S/H = 18.6 ppm. With the exception of S/H, 
these abundances were used by Mathis (2000) in his pho- 
toionization models of the local DIG. We found that the 
solar S/H (Anders & Grevesse 1989) provides a good match 
to the observations. The lower value, S/H = 13 ppm, used by 
Mathis (2000) was probably due to different dielectronic re- 
combination rates in his photoionization code (see discussion 
in § 3.5.4). We also investigate the effects of different abun- 
dances and perform some simulations using abundances ap- 
propriate for the Perseus arm (Mathis 2000) . We do not con- 
sider cooling from collisionally excited lines from elements 
other than those listed here, so our temperatures may be 
slightly hotter than ID models that include more elements 
(e.g., Sembach et al. 2000). 

Currently dielectronic recombination rates have not 
been calculated for third and fourth row elements, so pho- 
toionization modeling of sulfur lines are subject to uncer- 
tainties due to unknown atomic data (see Ali et al. 1991). 
We use the total recombination rates for S"*" and S^"*" from 
Nahar (2000). For S'^^ we use the radiative recombination 
rates from Verner & Ferland (1996) and the average dielec- 
tronic rate of 2.5 x 10~^^ from Ali et al. (1991). We discuss 
the issues related to modeling [S II] lines in § 3.5.4. 
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2.2 ISM Density Structure 

The ISM is observed to comprise several components of 
(warm and cold) neutral and (warm and hot) ionized gas. 
Unlike many recent photoionization models of the DIG, 
Monte Carlo techniques are not restricted to ID averaged 
models. For our smooth density models, we use the two com- 
ponent density from Miller & Cox (1993) 



n(H) = 0.1 exp(-|2|/0.3) +0.025 exp(-|2|/0.9) 



(1) 



where the number densities are per cm and the distances 
are in kpc. This represents the concentrated neutral layer 
and the extended ionized layer. In their model for the DIG, 
Miller & Cox used the above smooth density and included an 
approximation for absorption by dense clouds using a model 
that reproduces the statistics of clouds in the ISM. Although 
this smooth density is lower than the average density in- 
ferred for and H"*", when implemented in their "standard 
cloud" model using the known ionizing sources in the so- 
lar neighbourhood, it reproduces the average emission and 
dispersion measures observed in the Milky Way. 

An obvious criticism of the Miller & Cox model is that 
the density they used is smaller than that inferred for the H I 
in the Galaxy (e.g.. Dickey & Lockman 1990). However, they 
did use the known distribution and ionizing luminosity of O 
stars in the Solar neighbourhood, and found that the above 
density allowed for the gas to be ionized at large z along 
with reproducing many of the observations of the DIG. In 
reality the ISM is clumped on a large range of sizescales, so 
the Miller & Cox density provides an estimate of the smooth 
component required to allow ionizing photons to penetrate 
to large \z\. Therefore, such a density is a good starting 
point for 3D models that incorporate smooth and clumpy 
components: the smooth component should be close to that 
used by Miller & Cox. 



2.3 Ionizing Sources and Spectra 

We consider single sources within multi-component ISM 
density distributions. The ionizing spectra are taken from 
the WM-basic model atmosphere library (Pauldrach et al. 
2001; Sternberg, Hoffmann, & Pauldrach 2004) which pro- 
vides model atmospheres and emergent spectra including 
the effects of non-LTE line-blanketing and stellar winds. We 
did not run the WM-basic models, but used spectra from 
the library generated by Sternberg et al. (2004) . The spectra 
used were for model atmospheres with solar abundances, and 
gravities and temperatures in the range 3.6 ^ log g ^ 4.0 
and 30000 K < T* < 50000 K. 
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Figure 1. Total projected intensity and line ratio maps for Ha, 
[N II]/Ha, [S II]/Ha, and [S II]/[N II]. Axes are labeled in kpc. 
The anti-correlation of [S II] and [N II] with Ha and the rise 
of [N II] /Hq, [S II] /Hq towards the edge of the ionized zone is 
evident. [S II] /[N II] is fairly constant within the ionized region, 
but rises towards the outer edges. 
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Figure 2. Two dimensional slice through the centre of the simu- 
lation grid showing the neutral H fraction and temperature. The 
increase in temperature towards the edge of the ionized region 
due to hardening of the radiation field is clearly seen. 



3 STANDARD MODEL 

Our standard single source model has the two component, 
smooth density structure of Eq. 1, an ionizing luminosity 
Q(H°) = 6 X lO^^s"'^, and an ionizing spectrum for a WM- 
basic model atmosphere with T^, = 40000 K and logg — 
3.75. The output of the photoionization simulation is the 
3D temperature and ionization structure for the elements 
that we track. Using the output temperature and ionization 
grids, we then calculate emissivity grids for various lines and 



form maps of the intensity and line ratios. The results of our 
standard simulation are shown in Figs. 1-5. 

In forming the emissivity grid we do not include cells 
that have H^/H > 0.25. Emissivities for cells in our simula- 
tions that have higher neutral fractions are unreliable due to 
the resolution of our grid. Also, because we do not include 
the dynamics of ionization fronts, emissivities from such cells 
are not correct since the physics of the ionization front is not 
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included in our simulations (see Wood et al. 2004 for further 
discussion) . 

Numerical noise in the Monte Carlo simulations is most 
easily seen in the figures showing cuts across the line ra- 
tio maps in Fig. 3 and the [N II]/[S II] cut in Fig. 5. The 
"jaggedness" of the lines in these figures can be taken as 
an estimate of the numerical noise in forming the various 
line ratios. Noise is also seen in the [S II] /Ha vs Ha and 
[S II]/Ha vs Ha scatter plots in Fig. 5, with some points de- 
viating from the clear anti-correlation. Although these dis- 
crepant points appear to give a large spread in an otherwise 
tight anti-correlation, they represent only a few percent of 
all points plotted and are generally towards the edge of the 
ionized zone where our photon statistics are poorest. 

Figure 1 shows maps of the Ha intensity (in Rayleighs), 
and the [N II]/Ha, [S II]/Ha, and [S II]/[N II] fine ratios. 
In calculating the Ha intensity we have assumed that all 
the emission from the simulation grid is at a distance of 
2.5 kpc, roughly the distance to the Perseus arm (Haffner 
et al. 1999). This figure immediately shows the anticipated 
increase in [N II] /Ha and [S II] /Ha with increasing distance 
from the central ionizing source. The increasing line ratios 
arise from the hardening of the radiation field and subse- 
quent increasing temperatures at large distances from the 
source. 

For this simulation [S II]/[N II] ^ 0.6, similar to that 
seen in the local DIG (Hafi^ner et al. 1999) , and rises rapidly 
at the outermost edges of the ionized region. The rapid rise 
is explained by the relative ionization states of the elements. 
Close to the edge of the ionized region, N is mostly N"*" and 
then transitions rapidly to N'^, while S is transitioning from 
S^^ to (there is no in our simulations). Therefore, 
S'''/S increases and N^/N decreases at the edge of the ion- 
ized zone leading to the large [S II]/[N II] line ratios at the 
outer edges of the line ratio map (see also Bland-Hawthorn 
et al. 1997; Fig. 9). The rapid increase in line ratios from the 
interface in our simulations is generally not observed in the 
DIG and may refiect a shortcoming of our pure photoion- 
ization models for this interface. 

Figure 2 shows a slice through the grid in the x-z plane 
showing the ionization of H and the 2D temperature struc- 
ture. The ionized volume is extended perpendicular to the 
plane due to the increased density in the midplane. This 
resembles the "Stromgren volumes " seen in other 2D pho- 
toionization simulations. The increasing temperature with 
distance from the source is clearly seen in this figure. 

Figure 3 shows line ratios as a function of z for vertical 
cuts at a; = and 0.4 kpc from the centre of the intensity 
maps in Fig. 1. The rise of [N II] /Ha and [S II] /Ha with 
increasing \z\ is apparent as discussed above. [O III]/Ha 
decreases with height above the plane, while [O I] /Ha and 
[O II] /Ha increase with increasing \z\. The He I/Ha ratio de- 
creases with \z\ because the radiation reaching these regions 
has had some of its He-ionizing photons absorbed at lower 
\z\. For low density gas at 10* K, He I/Ha = 0.5He+/H+ 
(Osterbrock 1989), so the He I/Ha ratio probes the helium 
ionization and in turn the ionizing spectrum of the DIG 
(Reynolds & Tufte 1995). Figure 4 shows the incident spec- 
trum and the spectrum that reaches \z\ — 2 kpc, clearly 
showing the hardening of the radiation field and suppres- 
sion of the He-ionizing photons. 

In Fig. 5 we show a vertical cut across the centre of the 
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Figure 3. Vertical cuts at a; = kpc (solid) and x 
(dots) showing the variation of line ratios with z. 
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Figure 4. hicidcnt source spectrum (upper) and spectrum reach- 
ing \z\ = 2 kpc (lower). Compared to the source spectrum, the 
spectrum at large |2:| is harder in the H-ionizing continuum and 
has had its He-ionizing photons suppressed 



Figure 5. Vertical cut at a; = kpc showing the variation of 
[S II]/[N II] with z. The scatter plots show the correlation of 
[S II]/Ha with [N II]/Ha and the anticorrelation of [S Ilj/Ho; and 
[N II]/Ha with Ha. In the [S II]/Ha — [N II]/Ha plot, the sohd 
lines mark the range of observations in the local DIG (Haifner 
et al. 1999) and the squares show results from NGC 891 (Rand 
1998). 
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[S II]/[N II] line ratio map. As in Fig. 1 we see [S II]/[N II] ~ 
0.6, similar to observations by Haffner et al. (1999). Figure 5 
also shows some scatter plots as another way of displaying 
the trends seen in the intensity and line ratio maps of Fig. 1. 
These scatter plots display the values for all pixels in the line 
ratio images of Fig. 1. Such scatter plots (e.g., Haffner et al. 
1999) show the correlation of [S II] /Ha with [N II] /Ha and 
the anticorrelation of [S II]/Ha and [N II]/Ha with Ha. The 
solid lines in the [S II] /Ha — [N II] /Ha scatter plot show 
the range of observations in the local DIG (Haffner et al. 
1999) and the squares are those for NGC 891 (Rand 1998). 
The increase in [S II] /[N II] towards the edges of the ionized 
volume in our simulations is evident in the change of slope 
in the [S II] /Ha vs [N II] /Ha scatter plot. Note that the 
dots are appearing to turn around at very large values of 
[S II] /Ha, where the [N II] /Ha ratio is decreasing due to 
the N"*" —^ N° transition. Our grid resolution does not allow 
us to follow the decrease of [N II] /Ha smoothly to zero. 
The change of slope in the [S II] /Ha [N II] /Ha scatter 
plot is not present in data from the local DIG, but may be 
present in the NGC 891 data. Again, either the DIG is fully 
ionized or our models do not correctly treat emission from 
the ionized/neutral interface. The variation of [S II] /Ha and 
[N II] /Ha against Ha shows an anticorrelation with the line 
ratios being largest where the Ha emission is weakest, in 
agreement with observations. 

We now show some results and discuss the effects of 
varying source luminosity, ionizing spectrum, inclusion of 
extra heating, and ionization by photons from leaky H ll 
regions. The intensity and line ratio maps are qualitatively 
similar, so our discussion mostly focuses on the intensity 
cuts and scatter plots for the simulations below. 



3.1 Varying Ionizing Luminosity 

In Fig. 6 the ionizing luminosity is varied in the range 
2 X 10"^ s"^ < Q(H°) ^ 8 X lO"*^ s-^ All other parameters 
are the same as presented in the previous section. Increasing 
or (iccrcasing the source luminosity yields larger and smaller 
ionized volumes. For the lowest ionizing luminosity, the sim- 
ulation is radiation bounded at \z\ =0.7 kpc and the vertical 
cut shows that the [O I]/Ha ratio becomes large (~ 0.1) at 
the edge of the ionized volume. This is due to the combi- 
nation of high temperatures and oxygen rapidly becoming 
neutral at the edge of the ionized volume (see discussion in 
§ 4.2). For lower luminosity sources the [S II] /Ha ratio is 
quite large at large \z\. This is because the ionized zone is 
radiation bounded and we are seeing the effects of the inter- 
face (S^"*" — > S"*") described above. Beyond the edge of the 
ionized zone [S II] /Ha is formally infinite because the Ha 
intensity is zero. For this simulation, the resolution of our 
grid is not sufficient to see the rapid docroaso in [N II] /Ha 
that occurs in the transition zone where N"*" — > N". We just 
see the rise of [N II]/Ha and [S II]/Ha with \z\ and beyond 
the end of the dotted lines [N 11] /Ha = 0. 



3.2 Varying Ionizing Spectra 

The effects of varying the ionizing spectrum are displayed 
in Figs. 7 and 8. Line ratios and scatter plots are shown for 
WM-basic model atmospheres with solar abundances, and 
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Figure 6. Vertical cuts at x = showing the effects of vary- 
ing ionizing luminosity on line ratios. Ionizing luminosities arc 
2 X 10*'' s-'^ (dots), 4 X lO'^^s-i (dashed), 6 x 10"''' s^^ (solid), 
and 8 X 10"''' s^^ (dot-dashed). Low luminosities lead to a smaller 
ionized volume so there is a larger relative contribution from the 
ionized to neutral transition zone, hence some line ratios are much 
larger towards large \z\. The lowest luminosity source is radiation 
bounded at jzj = 0.7 kpc. 




Figure 7. Effect of varying ionizing spectra on line ratios ver- 
tical cuts at a; = kpc. Effective temperatures for the ionizing 

source are 35000 K (dots), 40000K (solid), 45000 K (dashed), and 
50000K (dot-dashed) . The harder ionizing spectra from the hotter 
stars lead to some very large line ratios. 



effective temperatures and loggr values of (35000 K, 3.8), 
(40000 K, 3.75), (45000 K, 3.9), and (50000 K, 4.0). The 

ionizing luminosity is Q(H'') = 6 x lO**^ s~^ and all other 
parameters are as in our standard model. As found in other 
investigations, harder spectra produce higher temperatures 
and increased line ratios for [N II]/Ha, [S II]/Ha, [O II]/Ha, 
and [O HI] /Ha. The 35000 K spectrum has little flux above 
24.6 oV in the Hc-ionizing contirmum, so the He I/Ha ratios 
are very small for this simulation. 



3.3 Composite Models and Extra Heating 

The results of Figs. 7 and 8 indicate that the observed 
line ratios in the local DIG and NGC 891 could be repro- 
duced with a range of source spectra and luminosities. We 
have followed Matliis (2000) and constructed models for a 
point source with a composite spectrum that has the fol- 
lowing contributions: 56% from T = 35000 K, 12% from 
T = 40000 K, and the rest from T = 45000 K. This is 
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Figure 8. Effect of varying ionizing spectra on scatter plots. 
Effective temperatures for the ionizing source arc (a) 35000 K, 
(b) 40000K, (c) 45000 K, and (d) 50000K. Tiie iiarder ionizing 
spectra from the hotter stars lead to increased line ratios, in some 
cases much larger than observed in the local DIG (solid lines) or 
NGC 891 (squares). 





Figure 10. Composite spectrum scatter plots (left) and with 
additional heating (right) Gi = 5x 10~^^ne ergs/s/cm^. The lines 
and squares show observations from the local DIG and NGC 891. 
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Figure 9. Composite spectrum models at a; = kpc (solid lines) 
and with additional heating (dots). 



comparable to the solar neighbourhood where the Garmany, 
Conti, & Chiosi (1982) catalog shows that around 50% of 
the O stars are spectral type 08 or later. Figure 9 shows 
that [N II] /Ha and [S II] /Ha increase with height above the 
plane, to values close to those observed in the local DIG. 
The line ratios are lower than those for NGC 891, likely 
reflecting the hotter sources ionizing its DIG. As with our 
other simulations, the very large [S II] /Ha line ratios from 
the ionized/neutral interface should be viewed with caution 
due to incomplete physics in our code. 

We have also investigated whether including additional 
heating may raise the [N 11] /Ha ratio for this composite 
model. Following Reynolds et al. (1999), we introduce an 
additional heating term, Gi, proportional to rie, which will 
dominate over photoionization heating at low Ue at large 
1 2: 1 in our simulations. The models in Fig 9 (dotted curves) 
and Fig. 10 have Gi — 5 x 10~^^ne ergs/s/cm^, yielding 
temperatures of around 11000 K at \z\ — 2 kpc. Increasing 
the heating to Gi = 10~^®ne ergs/s/cm"^, yields tempera- 
tures in excess of 15000 K at \z\ =2 kpc. The additional 
heating present in Figs. 9 and 10 is much lower than values 
inferred by Reynolds et al. (1999), Gi ~ 0.65-1.6 x lO'^'^ne 
ergs/s/cm^ and Mathis (2000), ~ 0.75-3 x 10"^^ ergs/s/H, 



Figure 11. Input spectra (a) and leaky spectra for leakage frac- 
tions of (b) 60%, (c) 30%, and (d) 15%. Note the hardening of the 
H-ionizing photons, the prominent 19.8 eV He I emission line, and 
the suppression of He-ionizing photons for small escape fractions. 

for ID models and is due to the natural increase of tempera- 
tures and line ratios at large |2| in our simulations. Figure 9 
shows the [S II], [N II], and [O II] emission is boosted com- 
pared to the simulation with no additional heating, with 
[O II] /Ha providing a very good probe of additional heat- 
ing (Mathis 2000). The very high [S II] /Ha line ratios from 
the ionized/neutral interface are still evident as discussed 
above. 

3.4 Leaky H ii Regions 

Observations indicate that many H ii regions are not tradi- 
tional radiation bounded Stromgren spheres, but axe density 
bounded and leak a fraction of the source luminosity into the 
DIG (Ferguson et al. 1996; Oey & Kennicutt 1997; Zurita et 
al. 2002). Depending on the source spectrum for an H ll re- 
gion, its He"*" zone may be smaller than the H"*" zone (e.g., 
Osterbrock 1989). Therefore, as suggested by Reynolds & 
Tufte (1995), if an H ii region is density bounded beyond the 
He"*" zone, but within the H"*" zone, the spectrum leaking into 
the ISM may have its He-ionizing photons suppressed and 
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Figure 12. Models showing the effects of Icalcy spectra on Une 
ratios vertical cuts at x = kpc. All models have the same ioniz- 
ing luminosity, but the ionizing spectra are as showrn in Fig. 11. 
The lines show a naked star (i.e., 100% leakage) 40000 K model 
(solid) and leaky models for escape fractions of 60% (dots), 30% 
(dashed) and 15% (dot-dashed). 




Figure 13. Effect of leaky spectra models on scatter plots. Panels 
show (a) naked star (i.e., 100% leakage), and models with escape 
fractions of (b) 60%, (c) 30%, and (d) 15%. 



the H-ionizing continuum hardened. This is seen in Fig. 4, 
where the spectrum leaking from a density bounded region 
is harder than the incident spectrum and may have its He- 
ionizing photons suppressed. Therefore, leaky H ii regions 
are a very plausible way for providing higher temperatures, 
due to the harder H-ionizing continuum, and lower He"*" frac- 
tions, due to the suppression of the He-ionizing photons. 
Hoopes & Walterbos (2003) presented ID leaky H ii region 
models, with some models producing elevated [S II]/Hq and 
[N H]/Ha and lower He I/Hq. We generate leaky spectra 
from a constant density H ll region with n(H) = 100 cm~^. 
The outer boundary is adjusted until the escape fraction is 
15%, 30%, and 60% of the incident source spectrum. The 
leaky spectra are then used as the source spectrum for 2D 
models of the DIG. All DIG models have the same luminos- 
ity, 6 X W*^ s~^, but differ in the source spectra, shown in 
Fig. 11. 

The incident (T* = 40000 K, logp = 3.75) and escaping 
spectra for our leaky H ll regions are shown in Fig. 11. The 
escaping spectra comprise photons that escape directly and 
contributions from the diffuse radiation field: H I Lyman 
continuum. He I Lyman continuum, He I two-photon con- 



tinuum, and the He I 19.8 eV, 2^5" ^ l^S emission hue. The 
leaky spectra are seen to be much harder in the H-ionizing 
continuum (flatter, or even increasing in the interval 13.6 eV 
- 24.6 eV) and for low leakage fractions the He-ionizing con- 
tinuum is suppressed. For leaky spectra with sources hotter 
than 45000 K, we found that the He-ionizing continuum is 
not significantly suppressed (see also Hoopes &i Walterbos, 
Fig. 13). 

Figure 12 shows the variation of line ratios with height 
for the leaky spectra of Fig. 11, and Fig. 13 shows the cor- 
responding scatter plots. The spectral shape of the leaky 
spectrum can lead to increased or decreased [S II] /Ha and 
[N II] /Ha (sec also Hoopes & Walterbos, Fig. 14). For the 
simulations shown here, [N II]/Ha increases as the H ll re- 
gion escape fraction decreases (all simulations have the same 
ionizing luminosity), while [S II] /Ha slightly decreases for 
the lower escape fractions. The increase of [N II] /Ha for the 
lower escape fractions is due to the harder spectra between 
13.6 eV and 24.6 eV for these models (Fig. 11). The He I/Ha 
decreases for leaky spectra due to absorption of He-ionizing 
photons in the H ll region. Thus some leaky H II regions pro- 
vide another way for producing the elevated temperatures, 
line ratios, and low He"*" fraction observed in the DIG. 



3.5 Other Models 

We have investigated many more models than those pre- 
sented above. Rather than present numerous figures, we now 
describe the main effects of varying other parameters. 



3.5.1 Varying Density Structure 

We performed simulations for a single component density 
structure, representing only the DIG component of the 
ISM. These simulations had lower temperatures, smaller 
[N II]/Ha, [S II]/Ha, and [O II]/Ha, and a larger He I/Ha 
than the two component model. This was due to the ab- 
sence of the concentrated density component which yields 
more hardening of the spectrum reaching large \z\ in the 
two component density models. Similarly, the plane parallel 
models presented by Rand (1998) and Bland-Hawthorn et 
al. (1997) did not produce as large [N II] /Ha and JS II] /Ha 
ratios as our two-component models. It appears that a uni- 
form slab, or single component exponential density struc- 
ture cannot yield the laxge [N II] /Ha, [S II] /Ha that are 
observed without additional heating of the DIG, or harden- 
ing of the radiation field via leaky H ll region models or a 
multi-component ISM density structure. 



3.5.2 Varying Source Location 

As the source location is moved to larger z heights above 
the plane, the ionized zone becomes radiation bounded for 
z < (e.g.. Miller & Cox 1993). At large positive z the grid 
is more ionized since the radiation has less of the concen- 
trated higher density component to ionize and pass through 
before reaching the gas at large z. Compared to a simulation 
where the source is in the midplanc, vertical intensity cuts 
for models where the source is located above the midplane 
show larger line ratios towards the lower radiation bounded 
edge and lower line ratios towards the upper edge of the 
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grid. These effects are seen in the simulations in § 3.1 where 
the source luminosity is varied. 

3.5.3 Varying Elemental Abundances 

Using abundances appropriate for the Perseus arm: (C, N, O, 
No, S)/H = (85, 45, 250, 71, 8) ppm (Mathis 2000) increases 
the [N II]/Ha, [S II]/Ha, [O II]/Ha, and [O m]/Ra line 
ratios by a factor of about 1.5 to 2 compared to those in 
Fig. 3. The lower abundances compared to our standard 
model result in increased line emission to obtain equilibrium 
temperatures, and hence the elevated line ratios. 



Some galaxies do show changes in slope of the [S II] /Ha 
vs [N II] /Ha and this has recently been interpreted by 
Elwert, Dettmar, & Tiillmann (2003) as an indicator for 
chemical evolution in galaxies. They suggest that increased 
[S II] /Ha compared to [N II] /Ha may arise from younger 
DIG layers since in the ISM nitrogen (from low mass stars, 
planetary nebulae, and stellar winds) is enriched more slowly 
than sulfur (from spuernovae type II). Our models adopt 
uniform abundances throughout the simulation grid and do 
not address this scenario. 



3.5.4 Sulfur Dielectronic Recombination 

Accurate modeling of the [S II] emission is not yet possible 
due to the unknown dielectronic recombination rates for S 
(see discussion in Ali et al. 1991). In our simulations we used 
total recombination rates for S"*" and S^^ from Nahar (2000) 
and the average suggested for S^+ from Ah et al. (1991). The 
rates of Nahar (2000) are from theoretical calculations, so 
may not be as accurate as from experiments. We have run 
our standard simulation with no dielectronic recombination 
for S and also using the averages suggested by Ali et al. 
(1991) for S+, S^+, and S^+. We find the [S II]/[N II] line 
ratio varies by about a factor of 1.5 among the three simular 
tions. Until more accurate dielectronic recombination rates 
are available, this uncertainty in modeling [S II] emission 
will remain. 



4 DISCUSSION 

4.1 Sulfur and Nitrogen 

The model line ratios presented above appear to be in good 
agreement with observations of [S II] and [N II] in the lo- 
cal DIG (Haffner et al. 1999), NGC 891 (Rand 1998), and 
several other galaxies (e.g., Otte et al. 2001, 2002). The in- 
crease of [N II] /Ha and [S II] /Ha with \z\ is a natural con- 
sequence of increasing temperatures away from the ionizing 
source due to the hardening of the radiation field. Compared 
to one dimensional averaged models, our two dimensional 
simulations reduce the requirement for extra heating to ex- 
plain the increasing temperatures and line ratios with height 
above the midplane. 

The change of slope in the [S II] /Ha vs [N II] /Ha scatter 
plots (e.g., Fig. 8) are not observed in the Milky Way's DIG 
(Haffner et al. 1999). The change of slope in our simulations 
is due to the increased S^/S and N'^/N fractions at the edge 
of the ionized volumes. The fact that such slope changes 
are not observed suggests that the DIG is almost fully ion- 
ized and not density bounded like our single source models. 
Wo will investigate multiple source models with overlapping 
ionized regions in a separate paper. Alternatively, it is quite 
likely that our models do not provide a good representation 
of the emission at the ionized/neutral interface as we do not 
include the effects of shocks or ionization fronts. The role of 
interfaces in interpreting ISM observations is very important 
(Reynolds 2004), and the very large [S II] /Ha and [N II] /Ha 
line ratios in our simulations may be a result of incomplete 
physics in our simulations. 



4.2 Oxygen 

There is currently limited data on oxygen lines in the Milky 
Way's DIG, with most observations probing the DIG at 6 = 
0°. However, in NGC 891 there are detailed observations 
of the dependence of [O III]/Ha and [O I]/Ha with height 
above the plane and Otte et al. (2001, 2002) have made 
measurements of [O II] /Ha in several galaxies. Rand (1998) 
finds that [O III] /Ha increases with height, which is opposite 
to what is seen in almost all of our simulations. The increase 
of [O III] /Ha with \z\ in NGC 891 is cited as evidence for 
O being ionized by a different mechanism, such as shocks 
(e.g., Collins & Rand 2001), instead of pure photoionization. 
Our models do not consider dynamics or ionization fronts, 
so cannot address these effects. Note, however, that models 
with very hard spectra (T* = 50000 K, Fig. 7) can produce 
increasing [O III] /Ha at large \z\. 

Otte et al. (2001, 2002) observe [O II] /Ha and 
[O III]/Ha to increase with height above the plane in five 
galaxies they studied, finding 0.5 < [O II] /Ha < 5. They 
also observed increases of [S II] /Ha and [N II] /Ha with \z\. 
Our pure photoionization models predict 0.5 < [O II]/Ha < 
3, and additional heating can raise this even further (e.g.. 
Fig. 9) . It appears that multi-dimensional pure photoioniza- 
tion models can reproduce most of the [O II] /Ha observa- 
tions of Otte et al. (2001, 2002), though additional heating 
or a harder ionizing spectrum may be required for some of 
the largest [O II] /Ha ratios. 

The [O Ij/Ha line ratios in our models show that [O I] 
increases in strength towards the edge of the ionized zone. 
The strong charge-exchange coupling of 0° to H° and the 
increased temperatures towards the edge of the ionized zone 
result in the increased [O Ij/Ha. This is generally seen as a 
problem with photoionization models, since in the few obser- 
vations in the local DIG, albeit at 6 = 0°, [O I] is observed 
to be rather weak with [O I]/Ha < 0.03. The [O I] emis- 
sion may be reduced if the region is fully ionized, or density 
bounded instead of radiation bounded (e.g., Mathis 2000; 
Sembach et al. 2000) . This is seen in the vertical cuts show- 
ing that [O I]/Ha < 0.03 in models where the gas is ion- 
ized beyond the maximum \z\ of our siniulatiou box (e.g.. 
Fig. 3). The increasing [O I] /Ha with \z\ in some of our 
simulations (e.g.. Fig. 6) do match the data for NGC 891, 
where [O I]/Ha ~ 0.1 at large \z\ (Rand 1998). More ob- 
servations of [O I] at larger \z\ in the Galactic DIG will 
determine whether there really is a difference in the [O I] 
emission between the Milky Way and NGC 891. 
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4.3 Helium 

As with [O I], tiie few measurements of He I in the Galactic 
DIG are at 6 = 0°. Tiie observations indicate that helium 
is underionized relative to hydrogen with He"'"/H+ < 0.04 
(Reynolds & Tufte 1995; see also Heiles et al. 1996). The He I 
observations probe the ionizing spectrum for the DIG and 
indicate a relatively soft spectrum, typically spectral type 
08 or later, corresponding to T* ~ 36000 K. The situation 
is not as extreme in NGC 891 where He"'"/H"'" ~ 0.06, indi- 
cating a harder ionizing spectrum for its DIG (Rand 1997). 

Our simulations using a composite spectrum, represen- 
tative of the O stars in the solar neighbourhood, produce 
fairly low values for He I/Ha, in line with current obser- 
vations at 6 = 0°. Further observations of He I at higher 
latitudes will provide additional tests of our models and the 
hardening of the radiation field at large \z\. 

We have also investigated an alternative mechanism for 
reducing the helium ionization, even when the sources are 
hotter than 08. The basic idea is that the radiation leak- 
ing out of midplane H ll regions to ionize the DIG may have 
its helium ionizing photons suppressed (e.g., Hoopes & Wal- 
terbos 2003). Our leaky H ii region models do indeed show 
that low values for He+ZH^ may be obtained, even if the 
ionizing source is quite hot. For the same ionizing luminos- 
ity, these models also produce higher temperatures than our 
other models, since the leaking ionizing spectrum is harder 
in the H-ionizing continuum. Therefore, leaky H II regions 
provide a plausible mechanism for explaining the low helium 
ionization in the Galactic DIG. 



5 SUMMARY 

We have presented some of the first multi-dimensional pho- 
toionization simulations investigating intensity and line ratio 
maps in models of the Milky Way's diffuse ionized gas. Our 
models produce [S II] /Ha and [N II] /Ha line ratios that are 
in good agreement with observations. Our simulations re- 
produce the observed increase of [S II] /Ha and [N II] /Ha 
with increasing \z\ due to the increasing temperatures at 
large distances from the ionizing sources. Previous analyses 
of line ratios in the DIG suggested additional heating was 
required to reproduce the elevated line ratios at large \z\. 
However, this appears to be due to the use of spherically 
averaged photoionization models. 

Most of our simulations do not include any additional 
heating above that provided by pure photoionization. When 
we did include extra heating it was much less than used 
in ID photoionization models of the DIG. Models such as 
those presented here, where line ratios are constructed for 
rays piercing the outer edges of photoionized regions, are 
more appropriate for the DIG. We have also shown that 
leaky H ll regions provide a plausible explanation for the low 
helium ionization in the DIG. The hardened leaky spectra 
produce higher temperatures than our standard models, in 
agreement with other work (Hoopes & Walterbos 2003). 

Unless we use very hard ionizing spectra, our models do 
not reproduce the observed increase of [O III] /Ha with \z\ 
in NGC 891. The observed rise of [O III]/Ha may be due to 
hotter ionizing sources or an additional source of ionization 
not present in our photoionization simulations (e.g., shocks 



or massive stars formed at largo \z\ from the dense clouds ob- 
served there, Howk & Savage 1997). Another potential prob- 
lem with our simulations is that they do not treat the physics 
at the ionization front, possibly resulting in an overpredic- 
tion of the [S II]/Ha line ratios. Despite those shortcomings, 
our models clearly demonstrate the importance of geometry 
in models of the DIG: compared to homogeneous or sin- 
gle component densities, multi-component models produce 
more hardening of the radiation field, and more elevated 
temperatures at large \z\. In a future paper we will explore 
3D models incorporating the known locations, luminosities, 
and spectral types of O stars in the Solar neighbourhood. 
The WHAM data in addition to the H I atlas (Hartmann & 
Burton 1997) should enable us to place constraints on the 
ISM density structure. 

We thank Alison Campbell, Torsten Elwert, Barbara 
Ercolano, Matt Hafi^ner, Kirk Korista, Lynn Matthews, John 
Raymond, Ron Reynolds, Jon Slavin, and Steve Tufte for 
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from an anonymous referee resulted in a clearer presentation 
of our results. KW acknowledges support from a PPARC 
Advanced Fellowship, JSM claims to be retired. 
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